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ABSTRACT 

The limited word size proposed for the GPS navigation computers could precipi- 
tate the problem of filter divergence in a standard sequential estimation algorithm. 
To insure filter reliability during the periods when the GPS observations are being 
processed, a square-root filter algorithm has been adopted for the GPS navigation 
computer. Several formulations of square-root filtering algorithms have been de- 
veloped during the past fifteen years. In implementing the filters, the covariance 
matrix can be propagated by either integrating the state transition matrix or by 
integrating a differential equation for the square root of the covariance matrix. 
The various approaches have different characteristics with re^rd to computer 
execution time, accuracy of the estimate, computer storage requirements, and 
the effort required to code and validate the algorithm. In this investigation the 
Potter, Carlson-Cholesky, and UDU square-root filters are compared with the 
standard extended Kalman filter. The characteristics of the algorithms are com- 
pared by simulating the application of a phase one GPS system to the determi- 
nation of a LANDSAT-D Satellite. 
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INTRODUCTION 

The limited word size used in contemporary microprocessor design 
may lead to problems in autonomous satellite navigation applications. 

The numerical error introduced when the navigation computations are 
performed with a short wordlength computer can lead to divergence of 
a standard extended sequential estimation algorithm. To insure filter 
reliability for applications where GPS observations are being processed, 
a square root filter algorithm has been adopted for the GPS navigation 
computer. Several formulations of square root filtering algorithms 
have been developed during the past fifteen years. This investigation 
describes a preliminary comparison of three square root filter formulations 
with the standard extended Kalman filter. Initial results are obtained 
regarding the relative computation speed and accuracy of these algorithms 
in a simulation of LANDSAT-D navigating with a Phase I GPS constellation. 

This summary is an overview of results obtained in the study. A 
complete discussion of the simulation procedure and numerical results can 
be found in [1] . 

t The research was supported by the National Aeronautics and Space 

Administration at Goddard Space Flight Center under Grant No. NSG 5154. 

1 

Professor, Department of Aerospace Engineering and Engineering Mechanics. 

2 

Research Assistant, Department of Aerospace Engineering and Engineering 
Mechanics. 

3 

Associate Professor, Department of Aerospace Engineering and Engineering 
Mechanics . 


254 



ALGORITHMS CONSIDERED 


The algorithms compared in this initial study are listed In Table 1. 
They are classified according to one of four ways of representing the state 
error covariance matrix and one of two methods of time propagation of this 
covariance matrix. 

The four methods of covariance representation are, briefly; 

1. standard Kalman formulation [2], 

2. the UDU algorithm [5], which decomposes the state error 
covariance into an upper unitary matrix, U, a diagonal 
matrix, D, and U'^^, 

3. the Carlson-Gholesky algorithm [4], in which the covariance 
is decomposed into an upper or lower triangular matrix, W, 
and its transpose an4 

4. the Potter algorithm [3] which decomposes the covariance 
into a general nxn square root covariance S and its trans- 
pose , s"^ . 

The two methods of perfoirming time updates of the covariance are 
referred to as the transition matrix method and the direct integration 
method. In the former technique, a system of variational equations 
is integrated from one measurement epoch to the next, to obtain a 
transition matrix. The time update of the covariance is performed by 
the proper multiplication of the covariance by the transition matrix 
and by the addition of process noise. The process noise matrix is 
obtained in this study by an approximate analytic integration of a 
diagonal spectral level density matrix. The direct integration methods 
involve the numerical integration of the state error covariance matrix 
directly. Direct integration is often considered for time propagation 
because, in general, fewer equations must be integrated between obser- 
vation epochs (n{n+l)/2 versus nxn for transition matrix methods). For 
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square root algorithms which use a transition matrix method, the in- 


clusion of process noise at each measurement epoch requires a retriangu- 
larization of the square root covariance matrix. This extra computati<»al 
burden is avoided in the direct integration methods as the process noise 
effects are included directly in the integration of the differential 
equations. 

The direct integration algorithm for the Potter filter requires 
an nxn matrix inversion at each integration step. As a result of the 
extreme numerical penalty resulting from the inversion operation, a 
directly integrated Potter algorithm is not considered. 


OBSERVATION SIMULATION PROCEDURE 

The filter algorithms were tested using a series of simulated range 
and range-rate observations made from the LANDSAT-D satellite to navi- 
gation satellites of a Phase I GPS constellation. The structure of the 
dynamic models used in the simulation of the observations is outlined in 
Table 2. The dynamic model of the motion of LANDSAT-D includes the 
geopotential effects of a non-spherical earth and the effects of 
atmospheric drag. The geopotential model used is GEM7 truncated to 
order and degree 8. The drag on the satellite is modeled as the function 
of a ballistic coefficient, the atmospheric density and the square of 
the user’s velocity relative to the atmosphere. The atmospheric density 
is determined from an exponential density model. 

In the observation simulation model, the GPS satellites are assumed 
to be in two-body circular orbits about a point mass earth. No other 
perturbations are assumed to affect them. 
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TABLE 1 


ALGORITHMS CtMPARED 

Transition Matrix Methods; 

• EKF(4) 

• UDU(I) 

• CARLSON-CHOLESKY(I) 

• POTTER (I) 


Direct Integration Methods: 

• EKF(P) 

• UDU(UD) 

• GARLSON-CHOLESKY(W) 

TABLE 2 

OBSERVATION SIMULATION MODEL 
User Model: 

• 8 X 8 Geopotential (GEM7) 

• Atmospheric Drag 

• Clock Errors 

GPS Satellites Model : 

• Two-Body Propagation Model 

• Clock Errors 
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The performance of the user's clock and the GPS satellites' clocks 
are important factors in determining potential navigation accuracy. The 
phase error for each of the clocks is modeled in the simulation as the 
sum of three terms: 

1. a deterministic phase error modeled as a first-order 
polynomial , 

2. the integral of an exponentially correlated frequency 
error, and 

3. a random walk phase error. 

The frequency error of each clock is the sum of derivatives of the first 
two terms of the phase error. 

Initial conditions for the Phase I GPS and for LANDSAT-D are shown 
in Table 3 and Table 4, respectively. The true anomaly of LANDSAT-D and 
the difference in the epoch times of the orbital elements can be varied 
to alter the pattern of GPS satellites visible from LANDSAT-D. The 
LANDSAT-D epoch elements are specified at a GPS system time of t = 0. 

The GPS epoch elements are specified at a GPS system time of -7200 seconds. 

Given these initial conditions, and given the dynamic models described 
previously, the history of GPS satellite visibility shown in Figure 1 was 
generated. Over the 21,500 sec. simulation period, the number of GPS 
satellites visible to LANDSAT-D varies from zero to six. From this 
visibility history, a set of simulated range and range-rate observations 
was determined for processing by the navigation simulation progr^mi. The 
observations were generated on the CDC6600/6400 system at the University of 


Texas at Austin. 



TABLE 3 


PHASE I GPS CONFIGURATION 


Satellite 

Long, of Asc. 
Node (Deg.) 

Mean 

Anomaly 

1 

-130. 

0. 

2 

-130. 

40. 

3 

-130. 

• 

o 

00 

4 

no. 

40 

5 

no. 

• 

o 

oo 

6 

no. 

120. 


Inclination: 63® 

Eccentricity; 0.0 

TABLE 4 

LANDSAT-D EPOCH ORBITAL ELEMENTS 

a = 7.086901 x lo^ m 
e = 0.001 
i = 98?181 
n = 354?878 
(0 = 180? 

f(true anomaly) = -185? 
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FIGURE 1 


GPS SATELLITE VISIBILITY VS TIME 



TIME (SEC) 
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FILTER STRUCTURE 


The liaite<l size of satellite on-board computer necessitates the 
use of simplified dynamic models in the navigation filter. Table 5 shous 
the dynamic model structure assumed by the filter for the user satellite 
and for the GPS satellites. 

The geopotential affecting the user satellite is of lower degree 
and order than that used to generate the simulated observations. The 
drag acceleration calculation uses an exponential model of atmospheric 
density, as does the observation generation program, but the ballistic 
coefficient of the user satellite is estimated. The bias of the user's 
clock at a given time is predicted by a linear equation. The coefficients, 
of this equation, the clock bias and drift at an epoch, are estimated. 

The estimation of the ballistic coefficient and of the clock parameters 
att«npts to account for modeling errors caused by the reduced size of 
the filter's geopotential. 

The filter prescribes a two-body point mass geopotential as the dyosmdc 
model for the GPS satellites. Each GPS satellite's clock bias is predicted 
between observation epochs by a linear equation whose coefficients are pre- 
determined to fit the clock's error behavior. For these simulations, the 
clock coefficients have been set to zero. This implies that the filter 
assumes perfect GPS clocks. Therefore, timing errors from all sources 
are accounted for in the estimation of the user's clock coefficients. 

Preliminary simulations have been performed with an eleven-state 
navigation filter. The differential equations of the states are shown in 
Table 6. The eleven states are: position (1-3) , velocity (4-6), ballistle 

coefficient (7), drag correlation parameter (8), user clock bias (9), user 
clock drift (10) , and clock drift correlation Parameter (11) . Important 
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TABLE 5 


FILTER SIMULATION MODEL 


User Model: 

» 4 X 4 Geopotential (GEM7) 

• Atmospheric Drag (Ballistic Coefficient Estimated) 

• Clock Error Polynomial (Coefficients Estimated) 

GPS Satellites Model: 

• Two-Body Propagation Model 

• Clock Errors (Predetermined Coefficients) 


TABLE 6 

FILTER DIFFEIffiNTIAL EQUATIONS 


11-STATE FILTER: 


r = V 


V 

d 



+ ^d 
d + 


^d 



(position) 

(velocity) 

(ballistic coefficient) 

(drag correlation parameter) 

(user clock bias) 

(user clock drift) 

(drift correlation parameter) 
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assumptions in the form of the differential equations are that the bal- 
listic coefficient and user clock drift are first-order Gauss-Markov 
processes and that the correlation parameters of the Markov processes 
are random walks. Additionally, the user clock's bias is modeled as the 
Integral of its drift. 

Parameters in the forcing functions of the equations in Table 6 are 
the following: 

a = gravitational acceleration 
a^ = drag acceleration. 

The are white noise forcing functions with statistics 

* 0 ; = Q. 

FILTER PERFORMANCE 

Figures 2 and 3 are plots of position error magnitude and velocity 
error magnitude of the estimates versus time. Two trends in the error 
history are evident from the examination of the plots. Catastrophic 
increases in the error magnitudes occur during periods of poor satellite 
visibility (less than four satellites in view) . A comparison of Figures 
1, 2. and 3 reveals the correlation between satellite visibility and 
error magnitude. The peaks of these errors reach over one kilometer 
in position, and over 1 meter per second in velocity. 

The second trend is observable by studying the long-term history 
of the errors. Although the errors decrease significantly when good 
satellite visibility is recovered, there is a long-term growth in the 
navigation error recorded during the periods of good visibility. As 
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an example, the position error before the first period of poor satellite 
visibility is at or below the 20 meter level. After the third period of 
poor visibility (approximately 16,000 sec.) position error has grown to 
approximately 100 meters. This indicates a secular or long period trend 
in the position error. A similar trend exists in the history of the velocity 
error magnitude. 

Although both error trends are significant and are deserving of 
further study, the goal of this Investigation has been to consider relative 
efficiencies of different algorithms. Therefore, current simulations have 
not been directed toward the removal of the error problems. This is to be 
the topic of follow-on studies. 


INITIAL NUMERICAL RESULTS 

The time propagation and measurement update algorithms described in 
the previous sections were tested on the 21,500 sec. arc of observations, 
using four different numerical integration schemes. The four schemes are; 

1) a variable step RK(7)8 with integration tolerances: 

^6 —10 

absolute error - 10 , relative error => 10 , 

2) a variable step RK(2)4 with integration tolerances; 

—2 —6 

absolute error = 10 , relative error » 10 , 

3) a fixed step RK8 with a step size of 6 seconds, and 

4) a fixed step RK4 with a step size of 6 seconds. 

Simulations with these integration algorithms are Intended to establish 
bounds on errors for navigation computations performed on the CDC6600 
system. The results of further studies using the lower order integrators 
proposed for the actual GPS computer can be compared to the results 
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shdva bere to determine the loss in navigation accuracy caused by using 
louez order Integration methods. 

iaStial navigation simulations, using die different filter algorithms 
over the entire 21,500 sec. observation arc, produced position error 
hlstofrles nearly identical to that In Figure '2. Specifically, no algorithm 
redvced the catastrophic position error growth during the third data drop- 
out Capproximately 12,000 sec.), or removed the long term error growth. 
Therefore, for the purpose of reducing computer use, the simulations for 
performance comparisons were run over the first 10,000 seconds of obser- 
vations only. 

The results of the initial set of simulations appear in Tables 7 
through 10. The pertinent values tabulated are: the total computation 
time for time update, the total computation time for measurement update, 
the total computation time (the sum of time update and measurement update 
cooptation times) , the total computation time normalized by the fastest 
total time, and the RMS error of the position estimate for the 10,000 
secoods of data. All computation times are tabulated in milliseconds, and 
RMS errors are written in meters. The algorithms are listed in the order 
of iiKreasing computation time. 

The results for the different filtering algorithms can be compared 
on the basis of computation time, estimation accuracy, and algorithm 
stability. Briefly, a comparison of the algorithms reveals the following 

trends. 

Computation Time 

The computation times of the algorithms display the fdllowing general 
character is t ic s : 
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TABLE 7 


RELATIVE ALGORITHM PERFORMANCE 
-6 -10 

RK(7)8 TOL: 10 - 10 


COMP. TIME (msec) 


ALG 

TIME 

UPDATE 

MEAS 

> UPDATE 

TOTAL 

UPDATE 

NORMED 

UPDATE 

EKF($) 

160.601 

9.435 

170.036 

1.000 

0DU($) 

172.119 

17.621 

189.740 

1.116 

POTT(I) 

179.924 

17.585 

197.509 

1.162 

CARL($) 

173.760 

24.968 

198.728 

1.169 

EKF(P) 

357.728 

9.442 

366.720 

2.157 

CARL(W) 

390.763 

24.772 

415.535 

2.444 

ODU(UD) 

530.734 

18.286 

549.020 

3.229 


ACC(m) 

iwi 

POS ERR 
128.4 
128.7 

95.0 

95.0 

129.0 
128.9 
128.4 



TABLE 8 

RELATIVE ALGORITHM PERFORMANGE 
RK(2)4 TOL: lO”^^ lO"^ 




COMP. TIME 

(msec) 


ALG 

TIME 

UPDATE 

MEAS 

UPDATE 

TOTAL 

UPDATE 

NORMED 

UPDATE 

EKF($) 

67.310 

10.215 

77.525 

1.000 

UDU(I) 

77.234 

17.629 

94.863 

1.224 

POTT($) 

85.569 

18.308 

103.894 

1.340 

CARL($) 

78.683 

25.256 

103.939 

1,341 

CARL(W) 

92.598 

24.560 

113.158 

1.550 

DDU(UD) 

144.782 

19.499 

164.281 

2.119 

EKF(f) 






ACC (m) 
RMS 

POS ERR 
128.4 
128.7 
95.0 
95.0 
128.3 
140.1 
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TABLE 9 

RELATIVE ALGORITHM PERFORMANCE 
RK 8, 6 SEC FIXED STEP 


COMP. TIME (msec) 


ALG 

TIME 

UPDATE 

MEAS 

UPDATE 

TOTAL 

UPDATE 

NORMED 

UPDATE 

EKF(I>) 

231.178 

8.957 

240.135 

1.000 

UDU($) 

243.286 

18.493 

261.779 

1.090 

POTT ( I) 

252.709 

17.675 

270.384 

1.126 

CARL(I) 

245.504 

26 .009 

271.513 

1.131 

CARL(W) 

444.624 

26.329 

470.953 

1.961 

UDU(UD) 

479.770 

18.191 

497.961 

2.074 

EKF(P) 

523.195 

8.801 

531.996 

2.215 


ACC (m) 
RMS 

POS ERR 
128.4 

128.7 

95.0 

95.0 

129.0 

128.8 

128.8 
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ALG 

EKF($) 

UDU(^) 

POTT(I) 

CASL($) 

CABL(H) 

EKF(P) 

UDU(UD) 


TABLE 10 


RELATIVE ALGORITHM PERFORMANCE 



M. 4, 6-SEC 

FIXED STEP. 




COMP. TIME (msec) 


ACC(m) 

TIME 

UPDATE 

MEAS 

UPDATE 

TOTAL 

UPDATE 

NORMED 

UPDATE 

RMS 

pos Err 

83.881 

8.408 

92,289 

1,000 

128,4 

95.613 

17.553 

113,146 

1,229 

128.7 

103.956 

17.884 

121.840 

1.320 

95.0 

96.959 

25.574 

122.533 

1.328 

95.0 

130.071 

25.632 

155.703 

1,687 

128.6 

153.496 

9.887 

163,383 

1.770 

128.9 



1) transition matrix algorithms consistently have lower computa- 
tion times for time update than do the direct integration 
methods , 

2) the relative performance of computation times of the direct 
integration algorithms is not predictable from one set of 
integration parameters to the next, contrary to the predicta- 
bility of the transition matrix methods, 

3) for this computer system, the transition matrix square root 
methods do not suffer large computation time penalties relative 
to the EKF($) algorithm, 

4) fixed step algorithms have larger computation times for time 
update and total computation time than do their variable-step 
counterparts. 

Accuracy 

Each algorithm generates an SMS position error of 128—129 meters, 
with tiK) exceptions; the UD simulation with the RK(2)4 integrator records 

an lafS error of 140 meters and the Carlson ($) and Potter ($) algorithms each 

• • 

record 95 meter SMS errors. The increased error of the UD simulation resultm 

—2 “6 

from the integration tolerances (10 and 10 ) not being strict enough to 

maintain the error at the 128 meter level. The integrator exceeds 
the error limit set for each step by the tolerances. The 95 meter SMS erroxs 
of the Carlson(l) and Potter(^) algorithms are thought to result from an 
assumption in the coding of the noise update sections of these two algorltlma* 
Currently, the reason for such a large drop in the position SMS is not knoen 
and, tterefore, the error improvement offered by these algorithms should be 
viewed with caution. 
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with the exception of the two cases Just mentioned, no algorithm 
pays a penalty in estimation error by decreasing the order of the 
numerical Integrator. 

Humerical Stability 

No algorithm appears to offer any clear cut advantage in stability 
in the problem conditions tested and for this computer system. However, 
the transition matrix equations do appear to be smoother and are, therefore, 
easier to integrate than are equations of the direct integration methods. 

This is evidenced by the generally lower integration times of the transitiam 
matrix codes and by the failure of two of the direct methods to successfully 
finish simulations. Both vailures resulted from the discovery of a negative 
valueii diagonal element in the covariance matrix brought about, it is bellsmed, 
by an inability of the integrator to maintain an adequate single-step error, 
thus causing a fatal global error growth. 


Additional Results 

In the previously discussed set of navigation simulations, the inte- 
gration of the state and state-error covariance matrix through data 
dropout periods is performed with one call to the numerical integrator. 
When only one integration call is made to span the entire data dropout, 
position error magnitudes during the dropout are not Included in the navi- 
gation program's calculation of the total RMS position error. As these 
errors should be included in the RMS calculation, the navigation program 
has been modified to integrate in six-second Intervals during dropouts 
and to calculate the position error at the end of each six-second interval. 
This means that numerical integration and process noise acciimulation are 
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performed in six-second increments until observations are obtained at 
the end of the dropout. As well as correcting the BMS position error 
calculation, this modification creates a better model of the operation of 
the actual GPS data processing system. 

A plot of the RSS position error, as generated by the modified 
navigation program, is shown in Figure A. Tables 11 and 12 contain the 
results of a set of navigation simulations in which the first 10,000 
seconds of observations only were processed. The simulations were per- 
formed with each of the algorithms tested previously and with the RK(2)A 
and RK4 integrators. The results given in Tables 11 and 12 are analogous 
to those contained in Tables 8 and 10, respectively, but were generated 
by the modified navigation program. 

A comparison of the data in Tables 11 and 12 with those of Tables 8 
and 10 shows no significant changes in the relative performance of the 
algorithms. The RMS position errors for most of the algorithms have grown 
from approximately 128 meters to approximately 166.5 meters because the 
position errors occurring during the data dropouts have been included in 
the calculation of the FmS error. The Carlson ($) and Potter ($) methods 
still have the lowest errors, but the values of the errors have experienced 
a growth roughly proportional to that experienced by the other algorithms. 
For the same reasons mentioned in the previous section, the lower RMS errors 
of these algorithms should be viewed with caution. Integration times 
and total computation times have Increased as well. The increases 
result because the newly imposed six-second integration interval in the 
data dropouts forces a step size limit on the variable step Integrator, 
causes more matrix retriangularizations for the transition matrix methods. 
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METERS -10 





TABLE 11 

RELATIVE ALGORITHM PERFOSMANCE 
RK(2)4, TGL: lO”^ - lO"® 



COMP. 

TIME (msec) 


ACC (m) 

ALG. 

TIME 

UPDATE 

MEAS 

UPDATE 

TOTAL 

UPDATE 

NOBMED 

UPDATE 

RMS 

POS ERR 

EKF($) 

79.296 

9.618 

88.914 

1.000 

166.5 

UDU($) 

91.916 

19.536 

111.452 

1.253 

166.5 

CARLSON ($) 

95.750 

26.607 

122.357 

1.376 

109.6 

POTTER ($) 

104.488 

19.361 

123.849 

1.393 

109.6 

EKF(P) 

119.353 

8.575 

127.928 

1.439 

166.9 

UDU(UD) 

225.471 

19.697 

244.968 

2.755 

166.9 

CARLSON (W) 

274.643 

26.735 

301.378 

3.390 

166.6 
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TABLE 12 


RELATIVE ALGORITHM PERFORMANCE 
RK4, 6-SEC FIXED STEP 



COMP. 

TIME (msec) 


ACC (m) 

ALG. 

TIME 

UPDATE 

MEAS 

UPDATE 

TOTAL 

UPDATE 

NORMED 

UPDATE 

RMS 

POS ERR 

EKF(I) 

91.196 

9.040 

100.236 

1.000 

166.5 

UDU($) 

104.366 

19.316 

123.682 

1.234 

166.5 

CARLSON ($) 

107,833 

26.065 

133.898 

1.336 

109.6 

POTTER (^) 

116,608 

19.320 

135,928 

1.356 

109.6 

CARLSON (W) 

136.829 

27.048 

163.319 

1.635 

166.5 

EKF(P) 

163.799 

8.652 

172.451 

1.720 

166.0 

UDU(UD) 

— 

... 

-- 


— 
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and generally creates a higher computation overhead due to a greater 
number of Integrator calls. 

The unexpected results in this new set of simulations are the Improved 

• • • 

performances of the P and UD algorithms in conjunction with the RK(2)4 
Integrator. In this set of simulations, both methods have RMS errors 
consistent with the errors of the other algorithms (approximately 166.5 

meters). Yet in the original simulations, those recorded in Table 8, 

• • • 

the P method failed to complete its simulation run and the UD code had a 

higher I^S error than did the other algorithms. The error Improvements 
are believed to occur because the step size limit' imposed by the six^second 

integration interval insures an accuracy in the integration of the P and 

• • 

UD equations that the integration tolerances cannot insure. 


SUMMARY 

It is reemphasized that these results apply to simulations on the 
CDC6600 system only. Although it is believed that these results do 
establish important trends in the relative performance of these algorithms, 
additional simulations on other comuter systems, particularly those with 
small wordlengths, are essential. Such simulations are currently being 
carried out. 

This article is intended as a general overview of the simulation 
method and initial results, A more exhaustive analysis of these results, 
as well as detailed descriptions of the simulation models and filter 
algorithms can be found in [1], 
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